#############################
#############################
###Spatial Mediation Model###
#############################
#############################
# Set RNG
set.seed(123)


#Prepare dat for analysis
# subset na rm data
med.dat <- na.omit(full.2003.2018[,c("gid","clear_zoon_confirm","lagclear_zoon_confirm","Forest.coverage_bin","ag_crops","fires.cell.count_anom","NDVI_mean_anom","log.nl","spllag_fires","spllag_ndvi","spllag_zoon","logged.pop","month","t")])
#Create regular dummies for the mediate package
unique_month <- unique(med.dat$month)
#run loop 
for (month in unique_month) {
  new_col_name <- paste0("f_month", month)  # Create a column name
  med.dat[[new_col_name]] <- ifelse(med.dat$month == month, 1, 0)  # Assign 1 if matched, else 0
}

#############
##Ag. areas##
#############
# subset only ag areas
med.dat.ag <- subset(med.dat, (ag_crops>0))

#Run base model
med.fit1.gs <- lm(spllag_ndvi ~ spllag_fires+spllag_zoon+
                 lagclear_zoon_confirm+log.nl +logged.pop+
                 f_month2+f_month3+f_month4+f_month5+f_month6+f_month7+
                 f_month8+f_month9+f_month10+f_month11+f_month12, data=med.dat.ag)
summary(med.fit1.gs)
#CSEs
med.cse1.gs <- coeftest(med.fit1.gs, vcov. = vcovCL(med.fit1.gs, cluster = med.dat.ag$gid, type = "HC0"))

#Run fit model
out.fit1.gs <- lm(clear_zoon_confirm ~ spllag_ndvi+spllag_fires+spllag_zoon+
                 lagclear_zoon_confirm+log.nl + logged.pop+
                 f_month2+f_month3+f_month4+f_month5+f_month6+f_month7+
                 f_month8+f_month9+f_month10+f_month11+f_month12, data=med.dat.ag)
summary(out.fit1.gs)
#CSEs
out.cse1.gs <- coeftest(out.fit1.gs, vcov. = vcovCL(out.fit1.gs, cluster = med.dat.ag$gid, type = "HC0"))

#Set seed
set.seed(123)
#Create a loop for simulations, accommodating memory limits
sim_results.gs <- list()
for (i in 1:10) {  # Split 1000 sims into 10 runs of 100 each
  sim_results.gs[[i]] <- mediate(med.fit1.gs, out.fit1.gs, 
                                 treat = "spllag_fires", 
                                 mediator = "spllag_ndvi",
                              sims = 100,   cluster = med.dat.ag$gid)
  gc()  # Garbage collection after each iteration
}

#Extract average p values
pvals_ag.gs <- med_pvals(sim_results.gs)

################
##Forest areas##
################
#Subset only forest data
med.dat.for <- subset(med.dat, (Forest.coverage_bin>0))

#Run base model
med.fit2.gs <- lm(spllag_ndvi ~ spllag_fires+spllag_zoon+
                 lagclear_zoon_confirm+log.nl + logged.pop+
                 f_month2+f_month3+f_month4+f_month5+f_month6+f_month7+
                 f_month8+f_month9+f_month10+f_month11+f_month12, data=med.dat.for)
summary(med.fit2.gs)
#CSEs
med.cse2.gs <- coeftest(med.fit2.gs, vcov. = vcovCL(med.fit2.gs, cluster = med.dat.for$gid, type = "HC0"))

#Run fit model
out.fit2.gs <- lm(clear_zoon_confirm ~ spllag_ndvi+spllag_fires+spllag_zoon+
                 lagclear_zoon_confirm+log.nl + logged.pop+
                 f_month2+f_month3+f_month4+f_month5+f_month6+f_month7+
                 f_month8+f_month9+f_month10+f_month11+f_month12, data=med.dat.for)
summary(out.fit2.gs)
#CSEs
out.cse2.gs <- coeftest(out.fit2.gs, vcov. = vcovCL(out.fit2.gs, cluster = med.dat.for$gid, type = "HC0"))

#Set seed
set.seed(123)
#Create a loop for simulations, accommodating memory limits
sim_results2.gs <- list()
for (i in 1:10) {  # Split 1000 sims into 10 runs of 100 each
  sim_results2.gs[[i]] <- mediate(med.fit2.gs, out.fit2.gs, 
                                  treat = "spllag_fires", 
                                  mediator = "spllag_ndvi",
                               sims = 100,   cluster = med.dat.for$gid)
  gc()  # Garbage collection after each iteration
}

#Extract average p values
pvals_for.gs <- med_pvals(sim_results2.gs)

######################
### All other areas###
######################
#Subset all other areas 
med.dat.other <- subset(med.dat, (Forest.coverage_bin==0|ag_crops==0))

#Run base model
med.fit3.gs <- lm(spllag_ndvi ~ spllag_fires+spllag_zoon+log.nl + logged.pop+
                 lagclear_zoon_confirm+
                 f_month2+f_month3+f_month4+f_month5+f_month6+f_month7+
                 f_month8+f_month9+f_month10+f_month11+f_month12, data=med.dat.other)
summary(med.fit3.gs)
#CSEs
med.cse3.gs <- coeftest(med.fit3.gs, vcov. = vcovCL(med.fit3.gs, cluster = med.dat.other$gid, type = "HC0"))

#Run fit model
out.fit3.gs <- lm(clear_zoon_confirm ~ spllag_ndvi+spllag_fires+spllag_zoon+
                 lagclear_zoon_confirm+log.nl + logged.pop+
                 f_month2+f_month3+f_month4+f_month5+f_month6+f_month7+
                 f_month8+f_month9+f_month10+f_month11+f_month12, data=med.dat.other)
summary(out.fit3.gs)
#CSEs
out.cse3.gs <- coeftest(out.fit3.gs, vcov. = vcovCL(out.fit3.gs, cluster = med.dat.other$gid, type = "HC0"))


#Set seed
set.seed(123)
#Create a loop for simulations, accommodating memory limits
sim_results3.gs <- list()
for (i in 1:20) {  # Split 1000 sims into 20 runs of 50 each
  sim_results3.gs[[i]] <- mediate(med.fit3.gs, out.fit3.gs, 
                                  treat = "spllag_fires", 
                                  mediator = "spllag_ndvi",
                               sims = 50,   cluster = med.dat.other$gid)
  gc()  # Garbage collection after each iteration
}

#Extract average p values
pvals_oth.gs <- med_pvals(sim_results3.gs)

#Combine all average p-values
pvals <- cbind(pvals_ag.gs, pvals_for.gs, pvals_oth.gs)
pvals

##Combine all estimates together
sim_results_list.gs <- list(sim_results.gs, sim_results2.gs, sim_results3.gs)

#Generate plots -- model name "main"
med_plots(model.name="Geospatial", sim_results_list=sim_results_list.gs)

##Export mediation models
# Open a file connection
sink("tableGS.txt")

# --- First: Export full regression tables (to LaTeX or text) ---
stargazer(med.cse1.gs, out.cse1.gs, med.cse2.gs, out.cse2.gs, med.cse3.gs, out.cse3.gs,
          type = "text", 
          star.cutoffs = c(0.05, 0.01), 
          style = "default")

cat("\n\n% --- Goodness-of-fit stats only ---\n\n")

# --- Second: Export only N, R², and Adjusted R² ---
stargazer(med.fit1.gs, out.fit1.gs, med.fit2.gs, out.fit2.gs, med.fit3.gs, out.fit3.gs,
          type = "text", 
          keep.stat = c("n", "rsq", "adj.rsq"),
          omit = ".*")

# Close the file
sink()


